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Abstract - The pairing properties of ultracold fermions, with an attractive interaction, loaded in 
a honeycomb (graphene-like) optical lattice are studied in a mean-field approach. We emphasize, 
in the presence of a harmonic trap, the unambiguous signatures of the linear dispersion relation 
of the band structure around half-filling (i.e. the massless Dirac fermions) in the local order 
parameter, in particular in the situations of either imbalance hoping parameters or imbalance 
populations. It can also be observed in the system response to external perturbation, for instance 
by measuring the pair destruction rate when modulating the optical lattice depth. Going beyond 
the mean-field level, we estimate the critical temperature for the "condensation" of the preformed 
pairs. 
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Since its first experimental observation, graphene has 
attracted considerable attention due to its interest in fun- 
damental physics as well as for potential applications [lp] . 
In particular, the energy band spectrum depicts "conical 
points" where the valence and conduction bands are con- 
nected, and the Fermi energy at half-filling is therefore 
only made of points. Around these points, the energy 
varies proportionally to the modulus of the wave-vector 
and the excitations (holes or particles) of the system are 
equivalent to ultra-relativistic (massless) Dirac fermions 
since their dispersion relation is linear [3] . In the presence 
of interaction, the vanishing density of states at the Fermi 
energy leads to extremely rich physics. Contrary to the 
square lattice, where the nesting of the Fermi surface gen- 
erally leads to ordered phases even for arbitrarily small 
interaction strengths, the Fermi-Hubbard model on the 
honeycomb lattice depicts two quantum phases transitions 
at half- filling 4 : first, for weak interaction, one observes 
a semi-metallic behavior with Dirac-like excitations J5j-[T] 
similar to the non-interacting situation; for an interaction 
strength U ~ 3.5J, the system enters into a spin-liquid 
phase, i.e. a Mott-Insulator (charge gap) without long- 
range magnetic ordering [8j|9j; eventually, for larger inter- 
action strength U ~ 4.3J, anti-ferromagnetic order sets 



Naturally, these very peculiar properties (massless Dirac 
fermions, spin-liquid...) have found a large echo in the field 
of ultracold atoms in optical potentials [lOj[Tl], which are 
now widely acknowledged to be great tools to study and 
understand the transport or magnetic properties of their 
condensed matter counterparts 12 13 . This is especially 



true in the case of graphene, since these quantum phase 
transitions have not been observed yet, due to possible 
discrepancies between graphene and the physics depicted 
by the Fermi-Hubbard model and to rather weak electron- 
electron interaction 12]; on the contrary the interaction 
strength between cold atoms can be tuned via Feshbach 
resonance. Therefore, even though atoms in cubic lattices 
is still the standard situation, cold gases in the honeycomb 
lattice start also to be the subject of many experiments, 



leading to new phenomena, either for bosons 14 15 or for 
fermions 
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Still, ultracold fermions in lattice have few drawbacks: 
the presence of an external confinement and a relatively 
high temperature 13 . Usually, the external confinement 



is taken into account using the local density approxima- 
tion, i.e. assuming that the local properties can be de- 
scribed by the non-confined system ones for the local 
chemical potential 17 -19 . On the other hand, this as- 



sumption usually breaks down when the fermionic density 
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has strong variation at the lattice scale; this is typically 
the situation at the border of the atomic cloud or, in the 
case of imbalanced population, at the boundaries between 
the different phases, like paired-non paired, polarized-non 
polarized... The proper description of these interfaces goes 
beyond the local density approximation. In this letter, the 
properties of a two-component fermionic gas with on-site 
attraction in the presence of both a honeycomb lattice 
and a harmonic trap are discussed in a mean-held ap- 
proach. To properly account for spatial inhomogeneities, 
the BCS order parameter at each site is an independent 
variable 
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whose value is determined, for a given 
temperature, by a global minimization of the free energy. 
After a short description of the method, we first describe 
the behavior of the order parameter in the standard sit- 
uation: perfect lattice and balanced population; then we 
emphasize the specificities of the honeycomb lattice when 
unbalancing either the hoping parameters or the popula- 
tions. Furthermore, we show how response functions of 
the system can be computed, for instance, the pair de- 
struction rate induced by a small modulation of the lat- 
tice potential. Finally, at large interaction, we emphasize 
the different temperature scales, namely the pair forma- 
tion and their condensation, the later transition being of 
the Berezinsky-Kosterlitz-Thoulcss (BKT) type, i.e. con- 
trolled by the phase fluctuations [22] . 

Starting from the full Fermi-Hubbard model on the hon- 
eycomb lattice (see Ref. 
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for a detailed description of 
the geometry) , one can derive the mean-field Hamiltonian: 



H MF = ^M^ + i J2 A* A.; - J2 



M 



-A* 
ft tt 



U 



(1) 



U{fLfJ,) are on-site pairing amplitudes; 



where A* 

* t = (■■■ ,fjp--- ,fil,--' J is the Nambu spinor. 
The matrix h depicts the one particle Hamiltonian, 
namely hopping terms between nearest neighbors hija = 
—Jij i 7^ j and chemical potential terms hu a = —fJ,i a = 
— Mct + ^Muj^pRf . In the case of perfect honeycomb 
lattice, all the tunneling rates have the same value 
J, which for a typical value Vq = IOEr, is given by 
J sa 0.07-Eh, where E R = h 2 k\/2M is the recoil en- 
ergy [II] . This value of J sets the unit of energy to 210 
nK (4.2 kHz) for 6 Li or 30 nK (600 Hz) for 40 Na. 

The free energy F = — \ In (Z) associated to the mean- 
field Hamiltonian reads: 



F 



^E A * A '-E^> ( 2 ) 



U 



where the are the 27V eigenvalues of the Nambu matrix 
M; N is the number of sites (2 x 51 x 51 in the following). 
The ratio between the harmonic trap energy fiwtrap and 
the bandwidth 6J is ps 0.04, i.e. corresponding to trap 
frequency 1.1 kHz for 6 Li, 160 Hz for 40 Na. 




Fig. 1: (Color online) Order parameter distribution in real 
space for two different values of the interaction U — 2.5J (a) 
and U = 5 J (b) (/3J = 25). In the first case, the order param- 
eter A depicts a dip having a ring-like structure corresponding 
to a local half-filling situation, at which, due to the conical in- 
tersections, the pairing is less efficient. At larger interaction, 
the ring disappear disappears and the A depicts a smooth pat- 
tern. The dip in the center reflects that, in general, the pairing 
efficiency decreases when getting closer to the unit filling. 



As explained in the introduction, the values of Aj are 
obtained by minimizing the free energy using a mixed 
quasi-Newton and conjugate gradient method; additional 
checks were performed to ensure that the global minimum 
has been reached. 

In the situation of balanced spin populations, the distri- 
bution of the order parameter is position space is shown 
in fig. [I] for two interaction values U = 2.5 J and U = 5J 
(/? J = 25). The value of the chemical potential at the cen- 
ter of the trap is fi = 2 J, leading to fermionic populations 
N a w 677 (resp. N a w 705) for U = 2.5 J (resp. U = 5J). 
For U = 2.5 J (a), the order parameter A depicts two dips, 
one in the center and one having a ring-like structure, cor- 
responding exactly to a local half-filling situation; this is 
nothing else but a signature that, at low interaction, the 
vanishing density of state at the the conical intersection 
leads to a poor pairing efficiency. At larger interaction, 
U = 5 J (b), this dip disappears and the A depicts a 
smooth pattern. The dip at the center arises since, for 
fermions in lattices, the pairing efficiency decreases close 
to integer filling. Therefore, when going from the outer 
part toward the center of the trap, the local number of 
pairs results from the competition between an increase of 
the number of fermions and a decrease of the pairing effi- 
ciency. Nevertheless, for a fixed trap geometry, the ratio 
between the total number of pairs, Yli |Ai| 2 /Z7 2 and the 
total number of fermion increases with the interaction: 
0.07 for U = 2.5J, 0.15 for U = 5 J, 0.24 for U = 10J and 
0.3 for U = 20 J. 



Imbalanced hopings. In Ref. 11 , it has been explained 
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Fig. 2: (Color online) Density and order parameter profiles 
along a diameter of the trap in the case of imbalanced hoping 
parameters (U — 5 J, /3J = 25). The black line with the circles 
corresponds to Aj (left axis). The red line with the squares 
corresponds to fermion density (right axis), (a): a situation 
where the gap is smaller than the interaction. The dip in the 
order parameter takes place around half-filling, where a kink 
in the density is already visible. For higher hoping imbalance 
(b), A now vanishes around half- filling, the density depicting 
a clear plateau. 



that the ratio between the tunneling rates can be con- 
trolled through the lasers building the optical lattice, al- 
lowing us to tune the band structure: the conical inter- 
sections can be moved, eventually merged, opening then 
a gap, i.e. leading to a two band model. This way, one 
achieves, at half-filling, the transition between the semi- 
metallic regime and a band insulator situation [9 , 23 



In 

the presence of an attracting interaction whose strength is 
smaller than the gap, one expects to observe the sequence 
superfluid-insulator-superfluid when increasing the chem- 
ical potential from below to above the band-gap. When 
the interaction strength is larger than the gap, the system 
will always be superfiuid, simply the order parameter is 
expected to become smaller when the chemical potential 
lies inside the band-gap. Indeed, in the situation when one 
of the hoping parameter J' is larger than the two others J, 
one observe, for J' = 2.5 J fig. [2^a), the presence of a dip 
in the order parameter around the position corresponding 
to half-filling, where the density depicts a kink, similarly 
to the non-interactive case 10 24 ; for higher hoping im- 



balance, it deepens, eventually leading, around J' — 3.25 J 
fig.^b), to a vanishing order parameter, i.e., splitting the 
superfiuid in two (independent) components. 

Imbalanced populations. The band structure of the hon- 
eycomb lattice is peculiar compared to the square lattice 
one: (i) there is no Fermi surface nesting at half-filling; 
(ii) the band structure is not separable, i.e. not given as 
the sum of two ID band structures. This results in a quite 
stable Sarma (breach pairing) phase, without a FFLO-like 
phase. For instance, for an interaction strength U = 5 J 
and a polarization P = 0.31, see fig. [3]ja), one recovers 



an usual situation 25 26 : the center of the trap is fully 
paired and not polarized, whereas the excess fermions are 
repelled on the edge, where one has a fully polarized phase. 
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Fig. 3: (Color online) Density and order parameter profiles 
along a diameter of the trap in the case of imbalanced popula- 
tion (U — 5 J, f3J = 25). The black line with the circles corre- 
sponds to A, (left axis). The red (resp. green) line with upper 
(resp. lower) triangles corresponds to the up spin (resp. down 
spin) density (right axis). The blue line with the squares cor- 
responds to the density difference, (a): polarization P = 0.31; 
the center of the trap is fully paired and not polarized; the 
excess fermions are repelled on the edge, where one has a fully 
polarized phase, (b): slightly higher polarization (P = 0.32) 
the order parameter now depicts a ring-like structure; the in- 
ner region of the trap is a partially polarized unpaired phase, 
the majority reaching a unit filling; the outer region is a still a 
fully polarized phase; the intermediate region is roughly made 
of a fully paired phase. 

For a slightly higher polarization P = 0.32, see fig. ^b), 
there is an abrupt change to the pairing distribution, with 
the formation of a ring-like structure of the pairs: the 
inner region of the trap consists in a partially polarized 
unpaired phase, the majority reaching a unit filling; the 
outer region is still a fully polarized phase; finally, the in- 
termediate region is roughly made of a fully paired phase, 
where the populations of the two fermion components are 
equal. This in a sharp contrast with the square lattice 
case, for which one would observe, for similar parame- 
ters, a checkerboard-like pattern for the order parameter, 
a FFLO-like situation 27 . 

Response functions. At the mean-field level, the 
Bogolioubov-De Gennes Hamiltonian being quadratic in 
the Nambu operators Vl/io-, all response functions, either 
for imaginary or real frequencies, can be computed from 
the single particle Green's functions in imaginary time, 
Giaj a i{j) = — ( v f/i cr (r)\l/] o .(0)), where i and j denote two 
lattice sites and (O) = Z- lr Tr[Oe- 0H ]. Denoting by P the 
matrix diagonalizing the Nambu matrix, i.e. M = P^AP, 
these single particle green's function have the following 
expression in imaginary frequency: 



G j. 



crjcr'K 



E 



iui n - A fe 



(3) 



where the cj n = 7r(2n + 1) denote the fermionic Mas- 
tubara's frequencies. 

For instance, one accessible experimental quantity is 
the pair destruction rate when modulating in time the 
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depth of the optical lattice, which translates into a time 
dependent modulation of the tunneling rate 



2N 



2\) 



Jij — 

Jj" + 8 cos (uit) . In the mean-field approach, the lo- 
cal average pair number is \^\ 2 /U 2 = (// t //j.) </ii/it) ^ 
such that the rate of pair creation/destruction is 
2Re((fyi)(f t (f ll f. n )}. The later quantity can be de- 
rived from the Heisenberg equations for the local pair op- 
erator MM'- 



-I 
'dt 



(MM)) = Moos (ojt) 4 ((MM) + (MM)) (4) 



the two average values (Mfjt) ano - (fjlM) a ^ time t can 
be obtained from the linear response theory, assuming that 
the modulation of the lattice is started at time t = 0. 



(MM)(t) 



5 

47T 

-icut 



e i(u)-u)')t _ ^ 

duJ '— 7 j^XilM u ') 

i{ui — a/) 

X4jt(- w ') 



i(u + a;') 

where Xiljfi m imaginary time, simply reads 



(5) 



XiWtW = E ^</*^)/itM&(°)/*r(o)) 
- V j° 



<p,q> 



<p,q> it 

G Plil(- T ) G jUl( T ) - G <zT4(~ t ) G jTpt( t ) 



(6) 



Inserting eq. Q in the preceding expression allows us to 
compute Xiijt^) as follows: 



J « w + ie + Afc _ Afc , ( n f(-^ k ) ~ 



E j ; 



P9 cj + ze + Afc - A fe < 



(n/(A fc ) - ^(AfcO) (7) 



where n/(A) is the Fermi function at temperature l/(/3J). 

From eq. Q and eq. ([5]), one easily sees that the pair 
rate, at long time t, has two components mainly oscil- 
lating at ±2w and one component which is almost time 
independent. In the limit t — > oo, the later is what would 
be given by the usual Fermi golden rule. Here, we com- 
pute the relevant quantities at a finite time t — 80/iJ -1 , 
larger than the typical timescale of the dynamics, but 
short enough not to resolve the discrete spectrum due to 
finite size of the system. The resulting pair destruction 
rate is displayed in fig. [4j right column, for two different 
site positions, whereas the left column depicts the single 
particle Green's function for a site at the center of the 
trap, i.e. lm(G il7 icr (uj)), for both U = (a) and U — 5J 
(b), with the same frequency resolution as for the pair 
destruction rates. Since at the center of the trap, the 
chemical potential is fj, = 2J, excitations, for U — 0, cor- 
responding to the conical intersections are located around 



lu = —2 J, which is indicated by the arrow. The peak se- 
ries at lower frequencies actually corresponds to the levels 
in the harmonic trap for the fermions having an effec- 
tive mass M* set by the curvature of the band at k = 0, 
namely M/M* = (2?r) 2 /9 x J/E R , leading to an effective 
harmonic trap frequency value wj* rap = u)t Tap yMjM* = 
Wtrap27r/3 x yj J/Er, which, in the present case, gives rise 
to effective harmonic levels separated by the energy 0.15J. 
Since the Green's function shown in the figure corresponds 
to a site almost at the center of the trap, only even lev- 
els can be excited resulting in a frequency separation be- 
tween the peaks equal to twice the effective harmonic trap 
frequency w 0.3J. Furthermore, one can check that the 
first peak precisely corresponds to the ground state whose 
energy is — 3 J + ?iu;j!j. ap , i.e. to a frequency transition 
ui = — /i — 3 J + tlki* rap = —4.85 J. Eventually, the peak 
structure disappears at higher frequencies, where devia- 
tions from the quadratic dispersion become more impor- 
tant. The situation is roughly similar in the interacting 
case, but with the addition of a clear BCS gap whose size 
is 2Ai. Its presence also shifts the excitations to lower 
energies; for instance the rough position of the conical 
intersection is now given by — y / /if~+~A| « —2. 6 J, as in- 
dicated by the vertical arrow. These features also show up 
on the pair excitation rates (fig |4jright column) : whether 
one considers a site at the center (c) or a site away from it 
(d), one can see the BCS gap, whose size is 4A 4 . For larger 
frequencies, one can also see the peaks corresponding to 
the harmonic levels. In between, the region with low exci- 
tation rates corresponds to the vanishing density of state 
around the the conical intersection; this last property is 
specific to the honeycomb lattice and obviously absent in 
the case of the square lattice. 

BKT transition. In the large interaction limit, there are 
two distinct energy scales: the pair formation (ksT « U) 
and the pair "condensation" (ksT ~ t 2 /U), which, in 
2D, is driven by the phase fluctuations of the pair order 
parameters A,, resulting in a BKT-like transition. 

Neglecting quantum fluctuations the partition function 
reads: 



V[A, A*]e"#^ l A °l 2 Tr L-^«HA,A*)1 (g) 



where A = (A 1; • • • , Ajv)- The saddle approximation in 
this integral leads to the usual BCS gap equations. In 
the strong interaction limit, the amplitude fluctuations are 
gaped, whereas the phase fluctuations are gapless (Gold- 
stone mode). Writing Aj = |Ai|e I(?ii , one fixes the ampli- 
tude |Aj| to its mean-field value and only accounts for the 



phase fluctuations 30 ■ 32 



,Tr 



(9) 



The BKT temperature can be estimated by expanding 
the preceding formula at the gaussian level [6j[l7}[T8}[33] : 



(10) 
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Fig. 4: (Color online) On the left, the single particle Green's 
function corresponding to a site at the center of the trap are 
plotted for U = (a) and U = 5 J (b). For U = 0, a dip 
in the excitation spectrum, due to the conical intersections, 
appears around cj — — 2J. For U = 5J, this dip moves to 
— \//uf + A? w —2. 6 J, because of the BCS gap Aj. The peaks 
structure corresponds to the harmonic levels of the fermions 
having an effective mass given by the curvature of the band 
around k = (see text). On the right, the pair destruction rate 
when modulating the lattice potential is plotted as a function 
of the modulation frequency. Plot (c) corresponds to a site a 
the center of the trap, whereas plot (d) corresponds to a site 
slightly away from the center. The BCS gap is clearly visible. 
More interesting is that the peaks due to the harmonic levels 
and the dip due to the vanishing density of states around the 
conical intersections are also visible. 



The coefficient 1^ describes the local spin stiffness and can 
be evaluated from the mean- field Green's functions Q. 
In the homogeneous situation, 1^ = /, the BKT critical 
temperature is the simply given by HbTbkt = f-f- In 
fig. [5] the left column shows the order parameter A along 
a diameter of the trap as a function of the inverse temper- 
ature P for two different values of the interaction U — 5 J 
(a) and U = 10 J (b). For a wide range of temperature, 
the profiles are essentially not modified: up to j3J = 5 
for U = 5J and j3J = 1 for U = 10J; In addition, as 
expected, for larger interaction, the order parameter van- 
ishes at higher temperature: f3J ~ 1 for U = 5J and 
(3J s» 0.5 for U = 10 J. Figure [5jc) shows the stiffness 
parameters for links along a diameter of the trap, for 
different values of the interaction U. As expected, it starts 
from a low value for U = 4 J, increases up to a maximum 
value around U = 7.5J and then starts decreasing again 
for larger value. Note that for U — 20J, even though 
the profile is different, the value is roughly the same as for 
U = 4 J. Finally, fig.^d) summarizes, as a function of the 
interaction U, the two temperature scales: the "classical" 
pairing Ta (red line with squares, left axis) and the "con- 
densation" (i.e. a quasi-long range order) of pairs Tbktj 
roughly evaluated from the spin-stiffness, (black line with 
circles, right axis). The blue dashed line is the "classical" 
pairing critical temperature obtained for the homogeneous 
system (no trap) having the trap center density. 
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Fig. 5: (Color online) The left column depicts the order param- 
eter along a a diameter of the trap for increasing temperature 
values (lower j3J values), for two different interaction strengths 
U = 5 J (a) and U = 10 J (b). For a wide range of tempera- 
ture, the profiles are essentially not modified: up to /3J = 5 
for U — 5J and j3J — 1 for U = I0J; Furthermore, the order 
parameter vanishes at higher temperature for larger interac- 
tion: 0J w 1 for U = 5J and /3J w 0.5 for U = 10J. (c): the 
stiffness parameters Iij for links along a diameter of the trap 
is plotted for different values of the interaction U. The maxi- 
mum value is reached around U = 7.5J, corresponding thus to 
the highest temperature of the pair "condensation", (d): the 
two temperature scales are plotted as a function of the inter- 
action U ; the "classical" pairing Ta (red line with squares, left 
axis) and the "condensation" of pairs Tbkt (black line with cir- 
cles, right axis) The blue dashed line is the "classical" pairing 
critical temperature obtained for the homogeneous system (no 
trap) with the density corresponding to the one at the center 
of trap. 



The preceding gaussian approximation does not account 
for the possibility of vortices in the phases, which are at 
the heart of the BKT transition. In order to cure this 
problem, one would have to compute the full partition 
function (|9j). However, it is numerically too expensive for 
large lattices. A good approximation consists in using 
the partition function of the effective XY model given by 



the 1^ distribution (10). We have checked that for small 
system, the agreement is fairly good. The signature of the 
transition can be measured in the momentum distribution 
of the pairs. For T < Tbkt, the quasi- long range order 
(AjAt) w |Af F ||Af F |(e 1 ^-^)) leads to a strong peak 
in the momentum distribution of the molecules, i.e. a 
quasi- "condensate" , see fig. [6] This peak broadens and 
disappear due to the vortex proliferation, for a value of 
the temperature j3J f=a 11 (U = 5J), in agreement with 
the critical value ^bTrkt ~ 0.09 J extracted from the spin 
stiffness, see fig. [5] From the experimental point of view, 
this "condensation" temperature corresponds to a ratio 
T/T F ss 0.02, a value T/T F w 0.2 is enough to the observe 
the pairing formation, thus, within experimental reach. 

In conclusion, we have presented pairing properties of 
ultracold fermions in a honeycomb lattice and in a har- 
monic trap with on-site attractive interaction. We have 
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Fig. 6: (Color online) Pair density in the quasi-momentum 
space for different values of the temperature /3, for U = 5J; 
(a): PJ = 13, (b): pj = 12, (c): f3J = 11, (d): pj = 10. The 
disappearance of the quasi-condensate peak around PJ = 11 
is in agreement with the critical value IcbTbkt ~ 0.09J ex- 
tracted from the typical spin stiffness value, see fig. [5] From 
the experimental point of view, this "condensation" tempera- 
ture corresponds to a ratio T/Tf ~ 0.02. 



shown that one of the main feature of the honeycomb lat- 
tice, namely the conical intersection, leads to unambiguous 
signatures in the pair distribution. For instance, at low in- 
teraction, a circular dip in the pairing distribution appears 
around the local half-filled situation. This dip can even be 
deepened by unbalancing the hoping parameters, eventu- 
ally leading to a two component superfiuid. In the case of 
imbalanced populations, the absence of Fermi surface nest- 
ing gives rise to ring-like structure of the pairing, rather 
than a checkerboard pattern (FFLO situation). Further- 
more, this vanishing density of state can also be observed 
in the spectrum of response functions like the pair creation 
rate when modulating the optical lattice depth. Finally, 
we have estimated the two type of transition temperatures, 
namely the pair formation and the pair condensation, em- 
phasizing that for actual experiments the observation of 
the paring is definitely within reach, whereas reaching the 
pair condensation regime still require a slight improvement 
in lowering the temperature. Among the possible exten- 
sions of the present work, one can mention the study of 



fermions in effective gauge fields 34 or the study of local- 
ization properties in the presence of disorder [2] . In addi- 
tion, genuine dynamics can be studied in a time-dependent 



Bogolioubov-De Gennes approach 35 
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